set.seed(1234)
library(splines)
library(ggplot2)
library(viridis)
library(patchwork)
library(hexbin)
library(reshape2)
library(ResistPhy)
library(ape)
library(posterior)
library(cmdstanr)
library(bayesplot)
library(latex2exp)

run_mcmc <- F

Ingest synthetic usage

We Restrict ourselves to a 20 year span Load usage data and simulation script

usage_df <- simulated_usage()
source("simu_phylo.R")

Set simulation parameters and simulate phylogenies

q_u <- 1.2
q_t <- -2.7

n_tip <- 200

gamma_sus <- 1/60.0*365.0

gamma_res_u <- gamma_sus + q_u
gamma_res_t <- gamma_sus + q_t

t0 <- min(usage_df$time)
tmax <- max(usage_df$time)
sim <- sample_phylo(gamma_sus, gamma_res_u, gamma_res_t, n_tip, 123456, 1234567) 

df <- sim$traj
phys <- sim$phys
mr <- sim$most_recent_samp

Visualise (t)

N_0 <- 3e6
incidence <- 1.5e5 / N_0
mu <- incidence+gamma_sus

time_grid <- seq(from=t0, to=tmax, by=1e-4)
t_change <- (t0 + floor(0.33*(tmax-t0)))

beta_func  <- function (t) if (t < t_change) 1.02*mu else mu*(0.98 + 0.06*(t-t_change)/(tmax-t_change))
N_func <- function(t) N_0*exp(log(2)*(t-t0)/(tmax-t0))
usg_func <- yearly_usg_stepfunc(usage_df$usage, usage_df$time)

beta_vals <- sapply(time_grid, beta_func)
beta_df <- data.frame(time=time_grid, beta=beta_vals)

usg_vals <- sapply(time_grid, usg_func)
usg_df <- data.frame(time=time_grid, usg=usg_vals)

N_vals <- sapply(time_grid, N_func)
N_df <- data.frame(time=time_grid, N=N_vals)

N_plt <- ggplot(N_df, aes(x=time, y=N))+
    geom_line() +
    xlim(t0,tmax) +
    labs(y=TeX(paste0("$","N(t)","$")))+
    theme_minimal()+
    theme(axis.title.x=element_blank(),
        axis.text.x=element_blank(),
        axis.ticks.x=element_blank())
beta_plt <- ggplot(beta_df, aes(x=time, y=beta))+
    geom_line()+
    xlim(t0,tmax) +
    labs(y=TeX(paste0("$","\\beta(t)","$")))+
    theme_minimal()+
    theme(axis.title.x=element_blank(),
        axis.text.x=element_blank(),
        axis.ticks.x=element_blank())
usg_plt <- ggplot(usg_df, aes(x=time, y=usg))+
    geom_line()+
    xlim(t0,tmax) +
    labs(x="Time",y=TeX(paste0("$","u(t)","$")))+
    theme_minimal()

f_plt <- N_plt/beta_plt/usg_plt
pdf("Figures/functions.pdf",4,6)
f_plt
dev.off()
#> png 
#>   2
f_plt

Visualise usage

usg_fun <- yearly_usg_stepfunc(usage_df$usage, usage_df$time)
ggplot(data.frame(time=time_grid, usage=sapply(time_grid, usg_fun)), aes(x=time, y=usage))+geom_line()+theme_minimal()

plot(ggplot(df, aes(x=t))+geom_line(aes(y=I_2),color="blue") +geom_line(aes(y=I_3),color="red"))

Preview the phylogenies run BNPR to see how that looks

plot(phys[[1]])
axisPhylo()

plot(phys[[2]])
axisPhylo()

Run MCMC

out <- infer_costs2(phys, 
                    mr,
                    usage_df$usage,
                    usage_df$time,
                    t0,
                    tmax,
                    1/60.0,
                    365,
                    n_iter=2000, 
                    n_warmup=2000,
                    model="nodecay",
                    K=60,
                    L=6.5,
                    gamma_log_sd=0.3,
                    stan_control=list(adapt_delta=.99,
                        max_treedepth=13,
                        parallel_chains=4,
                        chains=4,
                        refresh=10
                    ))                    
saveRDS(out, "~/simu_rbf.rds")
out <- readRDS("~/simu_rbf.rds")
print(out$converged)
#> [1] TRUE

Inferred dynamics & Rt

source("Figures/figure1.R")
pdf("Figures/figure1.pdf",6,6)
plot(fig1)
dev.off()
#> png 
#>   2
fig1

Pairs Plot

source("Figures/figure2.R")
#> Warning: The `size` argument of `element_line()` is deprecated as of ggplot2 3.4.0.
#> ℹ Please use the `linewidth` argument instead.
#> This warning is displayed once every 8 hours.
#> Call `lifecycle::last_lifecycle_warnings()` to see where this warning was generated.
#> Warning: `aes_()` was deprecated in ggplot2 3.0.0.
#> ℹ Please use tidy evaluation idioms with `aes()`
#> This warning is displayed once every 8 hours.
#> Call `lifecycle::last_lifecycle_warnings()` to see where this warning was generated.
pdf("Figures/figure2.pdf",6,6)
plot(fig2)
#> Warning: The dot-dot notation (`..density..`) was deprecated in ggplot2 3.4.0.
#> ℹ Please use `after_stat(density)` instead.
#> This warning is displayed once every 8 hours.
#> Call `lifecycle::last_lifecycle_warnings()` to see where this warning was generated.
dev.off()
#> png 
#>   2
fig2

Traces & ppcheck

s1 <- plot_traces(out)
pdf("Figures/s1.pdf",10,10)
plot(s1)
dev.off()
#> png 
#>   2
s1

s2a <- plot_ppcheck_At(out,1)
s2b <- plot_ppcheck_At(out,2)
pdf("Figures/s2.pdf",6,6)
plot(s2a)
plot(s2b)
dev.off()
#> png 
#>   2
s2a

s2b

Hyppeparameter pairs

s3 <- plot_hyperpar_pairs(out)
pdf("Figures/s3.pdf",6,6)
plot(s3)
dev.off()
#> png 
#>   2
s3